************************************************
***Title: robustness_randomized_inference.do
***Creators: Joelle Abramowitz, Shooshan Danagoulian, and Owen Fleming*
***Notes: This file produces the randomized inference figures. The first figure shows the estimate from each randomized iteration alongside the baseline estimate, and the second figure shows a density plot of the distribution of the randomized iterations.

*For questions, contact
*Owen Fleming
*hg3490@wayne.edu
************************************************


**********SETUP
clear all
set more off
set matsize 11000
set seed 8675309
use data/for_analysis.dta, clear


**********PRODUCE ESTIMATES
local top = 500
local yvar count
local tvar pollen_q4_ls

mat storage=J(`top', 2, .)

ppmlhdfe count pollen_q2_ls pollen_q3_ls `tvar' $weather, absorb(county year_month month_day) cluster(county) tolerance(1e-06)

mat storage[1, 1] = 1
mat storage[1, 2] = _b[`tvar']
 
foreach i of numlist 2/`top'  {
               
di _cont "."
               
shufflevar `tvar'

qui ppmlhdfe count pollen_q2_ls pollen_q3_ls `tvar'_shuffled $weather, absorb(county year_month month_day) tolerance(1e-06)
		
mat storage[`i', 1] = `i'
mat storage[`i', 2] = _b[`tvar'_shuffled]
               
drop `tvar'_shuffled
               
}

svmat storage
keep storage*
rename storage1 trial
rename storage2 beta

save data/random_inf, replace


**********PRODUCE VISUALIZATION AND EXPORT
*Iteration figure
use data/random_inf, clear
twoway (scatter beta trial if trial > 1, ///
                mlcolor("118 118 118") mlwidth(medium) msize(vsmall) mfcolor(white) ///
                graphregion(color(white)) legend(off) ///
                xtitle("Iteration") ///
                ytitle("Treatment effect") ylabel(-0.07(0.01)0.07)) ///
                (scatter beta trial if trial == 1, msize(medlarge) msymbol(D)), yline(0, lpattern(solid))
graph export results/robustness_randomized_inference_iteration.png, replace
			
*Density
use data/random_inf, clear
scalar baseline = beta[1]
scalar baseline_text = baseline+.014

kdensity beta, nograph generate(x fx)

line fx x, sort xline(`=baseline', lpattern(solid)) xtitle("Treatment effect") leg(off) ytitle("Density") text(15 `=baseline_text' "Baseline") lcolor("118 118 118")
graph export results/robustness_randomized_inference_density.png, replace


**********ERASE
erase data/random_inf.dta
